Quasiharmonic (QH) density functional lattice dynamics is routine nowadays to compute vibrational response of a crystal at finite temperatures.^[1,2] For instance, density functional perturbation theory (DFPT) gives accurate results for non-metallic substances,^[2–5] when the effect of explicit temperature-dependent phonon frequency and electronic contributions is less important. On the other hand, for metals, electronic excitation is added to lattice-ion motional part to find Helmholtz free energy. However, results based on quasiharmonic approximation (QHA) are accurate to a temperature T ≲ 0.8T_M, where T_M represents normal melting temperature of a substance under consideration. Close to the melting temperature, anharmonic effects, such as temperature-dependent vibrational frequency, vacancy contribution, and so on, are non-negligible and important for correct assessment of thermal properties of materials.

Although anharmonicity is directly included in first principles molecular dynamics (FPMD) studies, the implementation is restricted as the number of atoms that can be managed on small computing systems is relatively small.^[6,7] In this regard, density functional theory (DFT) can efficiently be used to overcome these difficulties and one can also sample complete Brillouin zone (BZ) over large q-vectors. This can be combined with different high-T models to find highly accurate results. For example, Wang et al.^[8–10] and later other researchers^[11–15] have employed free-volume theory to compute vibrational free energy. The essence of the theory is to replace true phonon density-of-state (p-dos) by a delta function centered at mean frequency of vibrating atoms. This classical mean-field potential (MFP) approach gives good results for high-temperature thermodynamics up to melting temperature, but does not account for the low-temperature quantum effect. Since the dominant high-T anharmonicity can be related to deviatoric thermal stresses, recently Carrier and Wentzocovitch^[16] and Wu et al.^[17] have shown that relaxing these thermal stresses can modify the equilibrium volume. This effective equilibrium volume can be deduced either by (i) minimizing self-consistent full QH phonon dynamics volume or (ii) by ensemble average of internal energy in order to account for renormalized frequency. Thus both these schemes fold-in the effect of anharmonicity into an equilibrium volume, and propose much improved QHA. The scheme produced excellent agreement for various thermodynamic properties for substances of geophysical importance.^[16,17]

The thermal perturbative approach is yet another popular approach to compute anharmonic free energy.^[18–20] To the lowest order in T, the anharmonic free energy for mono-atomic system can be written as , with volume-dependent parameter a(V) to be determined and the Boltzmann constant k_B. Wallace^[21] has shown that in the first approximation, the intrinsic anharmonicity can be included through explicit temperature-dependent phonon frequency ω(V,T). Based on this assumption, Oganov and Dorogokupets^[18] have derived analytical expressions for various thermodynamic variables for weakly anharmonic oscillator. They have tested and validated the proposed scheme (within the Einstein model) to MgO combining to MD simulations. Based on this study, we infer that the use of the method is computationally very favorable up to melting temperature, and hence offers a very reliable scheme. In this paper, we study the role of anharmonicity through the lowest order for anharmonic free energy in conjunction with full QH-lattice dynamics. However, the present approach differs from the one given in Ref. [18] on two fronts. First, in order to derive analytical expressions for various thermodynamic properties, Oganov et al.^[18] have used Einstein model for phonon dynamics whereas we retain DFT-based full lattice dynamics. Second, though non-metallic MgO is prone to anharmonicity, the role of electronic excitations (due to strong sp–d hybridization) in d-band metals (such as Ag in the present case) makes it more sensitive to the external influence like temperature. We have combined first principles DFPT results to quadratic-temperature-dependent perturbation term for inclusion of anharmonicity. Parameters of this term are optimized through thermal expansion coefficient. We have calculated various thermo-physical properties up to melting temperature, taking silver as a test case.

The main objective of the paper is to report the thermodynamic results from a systematic inclusion of phonon anharmonicity in fcc-Ag. The next section gives overview of the computational methodology. We present comparisons of first principles quasiharmonic results and those derived after including anharmonic term to infer the role of anharmonicity at elevated temperature in Section 3. Based on the computed thermophysical properties and their comparisons with experimental and other theoretical results, we confirm the importance of the scheme in determining anharmonicity of higher order at high temperatures in the final section.

The Helmholtz free energy is given as a sum of cold energy E_C(V), QH lattice-vibrational energy , electronic free energy F_el(V,T) and anharmonic energy F_An(V,T), that is 1 Binding or cohesive energy E_C(V) is calculated as follows. The ab initio calculation for phonon dynamics were performed using PWSCF and PHONON code^[22] to derive QH free energy. The ion–valence electron interactions were described using an ultrasoft scalar-relativistic pseudopotential. The used potential also includes non-linear core correction, which was generated using Vanderbilt code with [Kr] 4d¹⁰ 5s¹ as the core state. A plane wave basis set with a cut-off of 120 Ry was used after checking the convergence for electronic SCF. To tackle the expansion of the augmentation charges, a relatively high cut-off of 800 Ry was used. We have used generalized gradient approximation (GGA) due to Perdew–Burke–Ernzenhof (PBE) for exchange-correlation functional. We prefer to use PBE as it is believed that the GGA gives better description of linear response of uniform electron gas compared to local density approximation (LDA), however improved exchange-correlation functional is also available.^[23] Total energies were then calculated self-consistently within DFT. Dynamical matrices were computed using 4 × 4 × 4 q-point mesh while dense 14 × 14 × 14 regular k-point mesh was implemented to Fourier interpolate these matrices. Force constants so deduced are used to compute p-dos and hence the other thermodynamic quantities. Since the present study aims to the metallic system, Methfessel–Paxton (MP) scheme with smearing parameter equal to 0.06 Ry was used. The total energy (E_tot) at various volumes were thus obtained by solving Kohn–Sham energy functional. Total energy at various volumes so deduced is then fitted to third-order Birch–Murngahan (BM) EOS to derive equilibrium volume V₀(≡ lattice constant, a₀), bulk modulus B₀, and its first order pressure derivative B′. These are shown in Table 1. These equilibrium properties are in excellent agreement with experimental and other theoretical findings.^[24–34] To derive cold energy curve, we need to calculate energy of an isolated atom. To achieve this end, we took large volume (V = 10V₀) such that practically it represents an isolated atom. Now using E_C(V) = E_tot(V) − E_atom, we have obtained cohesive energy at different volumes as depicted in Fig. 1.

Table 1.

Table 1.

Table 1. Equilibrium properties for fcc Ag. . Present Experiment Others a0/a.u. 7.6887 7.7289[24], 7.7365[25] 7.5967, 7.8669, 7.6780, 7.8518, 7.6840, 7.8934, 7.9728, 8.0540, 8.0805, 7.7550, 7.8159, 7.8272, 7.5117, 7.7252, 7.752 – all quoted from Ref. [25] B0/GPa 107.4 100.7[24], 108.0[26] 126.2[21], 117.4[27], 138.8[28], 117.0[29], 71.0[30], 102.0[29] B′ 4.251 4.11[31], 6.18[32] 4.725[33], 4.846[33], 4.270[34], 3.650[34], 4.540[34], 3.720[29]

Table 1. Equilibrium properties for fcc Ag. .

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) Cohesive energy curve. Present DFT results (solid squares) and BM-EOS fitted results (line) are compared with experimental result (star) due to Kittel.[24]

The lattice vibration part of free energy according to quasiharmonic approximation is given by 2 where wave-vector q and subscript j index the lattice vibration frequencies ω_q,j(V). The first term represents zero-point energy. In practice, equation (2) is conveniently converted to integration over the first BZ through the concept of p-dos, normalized to number of atoms. In the present study, phonon dispersion curve (pdc) and associated p-dos are calculated with the DFPT. Figure 2 shows the pdc and p-dos for equilibrium volume. Result for pdc is in good agreement with experimental results^[35] at T = 0 K. We have computed pdc and p-dos for several lattices constant ranging from 6.50 a.u. to 8.50 a.u. in the step of 0.05. The QHA module implemented in Quantum Espresso code^[22] is then used to deduce at various volumes and for several temperatures. This is needed to minimize free energy in determining equilibrium volume V₀(T) at each temperature.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Phonon dispersion curve and phonon-density-of-state. Experimental results (solid square) are from Ref. [35].

The electronic part of the free energy is given by 3 where E_el and represent internal energy and entropy, respectively. In the present study, E_el(V,T) is evaluated using semi-empirical formula due to McCloskey^[36] 4 where , is the coefficient of electronic specific heat, and with gas constant R, atomic number Z, and . For silver, Z = 47. The coefficient of electronic specific heat depends on the volume through the relation . Here, erg⋅gm⁻¹⋅K⁻² taken from Ref. [36]. It is proposed^[36] that the equation (4) can correctly interpolate high temperature ideal electron gas limit to the low temperature degenerate electron gas formula, viz. .

As discussed in the Introduction, F_An(V,T) is written in terms of lowest order term in temperature 5 According to Refs. [18]–[21], the physical meaning of parameter a(V) can be Since it is impossible to compute explicit temperature dependence of phonon frequency we treat parameter a(V) as a fitting parameter. Its volume dependence^[18] is usually given by 6 where constants and m are known as anharmonicity fitting parameters. Since the principal effect of anharmonicity is on the thermal expansion mechanism, we have fitted parameters and m to linear thermal expansion coefficient close to melting temperature, where one expects largest anharmonicity. We choose T = 1000 K to optimize these parameters. The melting temperature of Ag is 1234.8 K. We can identify the sum as the total (including anharmonic contribution) vibrational (≡ lattice) free energy.

These four contributions are added to compute total Helmholtz free energy in Eq. (1). Equilibrium volume V₀(T) and corresponding lattice constant a₀(T) at each temperature are obtained by minimizing F(V,T). In achieving this end, we have used 29 volume points ranging from V = 0.7V₀ to V = 1.4V₀ in the step of 0.025 to calculate total energy and hence the cohesive energy E_C(V). This is shown in Fig. 1. Also, we have performed DFPT calculations to deduce pdc and p-dos for each of these 29 volumes, and hence QH lattice vibrational free energy at each temperature (0–1300 K) using QHA module implemented in QE-code.^[22] Similarly, the electronic free energy in Eq. (3) is also calculated for the same volumes and temperatures. The anharmonic contributions, Eq. (5), are evaluated with some guess value of and m with Eq. (6) at each temperature and volume. This tabular information is then added to give the total free energy versus volume for each temperature. Fitting third order BM-EOS to each energy–volume curve gives equilibrium volume V₀(T) at an assumed temperature. Based on these V₀(T), we have first computed linear thermal expansion coefficient We then tune parameters m and and repeat this computational cycle till we get the best fit for α_L (T = 1000 K) to experimental results. Finally, this gives and . These refined values of and m are used to recalculate F_An(V,T) and hence the total free energy at each temperature. F(V,T) is then again fitted to BM-EOS to yield final equilibrium V₀(T), isothermal bulk modulus B(T), and its first pressure derivative B′(T).

Since B_T is now available through BM-EOS fit to total free energy at each temperature, adiabatic bulk modulus is obtained as 7 where C_V denotes constant-volume specific heat, which is discussed below. Total C_V can be written as the sum of the lattice-ion, anharmonic, and electronic contributions, that is 8 where 9 10

The electronic and anharmonic contributions have been computed by Eqs. (4)–(6). Lattice contribution is available from QHA module.^[22] Furthermore, the thermodynamic definition of Grüneisen parameter is Now, theoretical constant pressure specific heat C_P is given by 11 According to Burakovsky and Preston,^[37] the thermodynamic Grüneisen parameter (γ_th) at ambient condition can also be calculated using 12 where B_T and denote isothermal bulk modulus and its pressure derivative, and P is pressure. The parameter λ(V) ≡ λ = −1.0, 0, 1, and 1.5 corresponds to expressions for γ_th due to Slater,^[38] Dugdale and MacDonald,^[39] Vaschenko and Zubarev,^[40] and improved free-volume theory by Stacy,^[41] respectively. Here, we assume that the equation (12) is also valid at different temperatures with corresponding equilibrium quantities. Moreover, it was also found in previous studies^[10] that the choice λ = −1.0 gives better thermal expansion coefficient. Thus, with presently estimated B_T and at each temperature under zero-pressure condition, we have calculated γ_th(V₀(T)) with λ = −1.0.

Present ab initio results for cold energy E_C(V)(=E_tot(V) − E_atom) as a function of volume are shown in Fig. 1 along with BM-EOS fitted results. We obtain cohesive energy E_C = −0.217 Ry in agreement with experimental results −0.216 Ry^[24] and −0.218 Ry.^[25] Computed equilibrium properties, such as lattice constant, bulk modulus, and its pressure derivative, are in excellent agreement with experimental and other findings as reported in previous section, shown in Table 1.

Since the present approach is based on first principles full lattice dynamics supplemented by empirical anharmonic term, we compare our phonon results in Fig. 2. One can see a very good agreement with experimental findings due to Kamatikahara and Brockhouse.^[35] Similarly, we have calculated pdc and p-dos (not shown here) for each volume to deduce ionic part of free energy at each temperature. This rigorous approach allows us systematic study of the quasiharmonic and anharmonic contributions to vibrational free energy.

Figure 3 represents lattice part of free energy versus temperature. Though both graphs appear similar, they differ in magnitude. However, as T → 0, anharmonic contribution will also go to zero. At T = 0 K, they converge to positive zero-point energy of 1.486 mRy. Furthermore, it is found that at temperatures 185 K for QHA and 165 K for total (including anharmonic contribution, written as An), the lattice free energy is zero, and crosses from positive to negative value. It is interesting to identify this temperature as the characteristic temperature (T_C) separating quantum and classical regime. These values are close to conventional Debye temperature for Ag, θ_D = 225 K.^[24]

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) Lattice vibrational free energy (blue: An, red: QHA) as a function of temperature.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) Linear thermal expansion coefficient. Present results: lines (blue: An, red: QHA). Experimental result (symbol) is from Ref. [43]. Theoretical results:[44,45] dotted line (QSC), dashed line (MSC), and dashed-dotted line (SC).

Moreover, using the corrected equilibrium volumes V₀(T) as described in the previous section, we have calculated the linear thermal expansion coefficient, shown in Fig. 4. Again, mutual comparison between QH-based results and the one obtained after correcting for anharmonic term reveals the importance of anharmonicity for better thermal expansion. Present anharmonic results are in excellent agreement with experimentally recommended values by Tolukin et al.^[43] Shown in Fig. 4 are the theoretical results due to Januszko and Bose.^[44,45] These authors have used different versions of the Sutton–Chen (SC) potential, viz. modified SC (MSC) and after including quantum correction to SC potential (QSC), to assess the applicability of these potentials at high-T. Mutual comparison of these data shows large scatter in thermal expansion coefficient, and indicates large differences in elastic strain produced by these potentials. Too high results actually reveal too-soft behavior compared to experimental trend and give unacceptable results.

Isothermal and adiabatic bulk moduli are depicted in Figs. 5 and 6, respectively, where comparison is also made with reported findings.^[44–47] We can infer that all results decrease with temperature. However, close to the melting temperature, bulk moduli decreases rather rapidly and non-linearly to account for rapid thermal expansion. Again, theoretically derived bulk moduli in Refs. [44] and [45] shows a large scatter in bulk moduli when obtained by three versions of SC potentials, suggesting different elasticity.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) Isothermal bulk modulus at different temperatures. Present result: lines (blue: An, red: QHA). Theoretical results: green dashed line (MSC),[44,45] black dashed line (SC),[44,45] purple long-dashed line (QSC),[44,45] red square with connecting line,[46] and green triangle with connecting line.[47] Experimental results (blue star) are from Ref. [21].

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. (color online) Adiabatic bulk modulus. Present result: lines (blue: An, red: QHA). Theoretical results:[46,47] magenta dashed line (MSC), solid black line (SC), and black dashed line (QSC). Experimental results (red down triangle) are from Ref. [21].

Figure 7 shows the results for C_V (T) within the QH (without term in Eq. (8)) and full anharmonic approximations, Eq. (8). Total phonon contribution ( remains an order of magnitude larger than the electronic ones (not shown here). The latter increases linearly with temperature, which is significant only close to melting temperature. Furthermore, harmonic approximation saturates to Dulong–Petit (DP) law (3k_B) at higher temperatures, whereas anharmonic contributions is near-linear in temperature with negative slope 3k_Ba(V). On this ground, according to Ref. [49], it is possible to classify the crystal into two types: the one for which is positive and the other one with negative value. Stacey et al.^[49] have shown that crystals with bonds to each atom appearing in opposite pairs, i.e. symmetrical, have negative values for . These crystal structures are simple and closed packed, e.g. fcc. In unison, Zoli^[50] has found that for aluminum and copper is negative. Zoli found that is −0.44 J⋅mol⁻¹⋅K⁻¹ for copper at T = 300 K, and similar (in magnitude) results can be expected to other noble metals as well. For example, we have also calculated for copper using the same computational scheme, and find that its value is −0.304 J⋅mol⁻¹⋅K⁻¹. The present estimate for anharmonic contribution is −0.74 J⋅mol⁻¹⋅K⁻¹ at T = 300 K. The same conclusion can also be drawn from previous findings.^[8,9,12] Thus decrease in C_V at higher temperature is the manifestation of anharmonicity. The larger the departure from DP law is, the stronger the anharmonicity is.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. (color online) Variation of specific heat at constant volume. Present result: lines (blue: An, red: QHA). Theoretical results:[44,45] red square (Morse potential) and black square (modified Morse potential). Horizontal dashed line indicates Dulong–Petit law.

The present results for C_P are compared with other findings in Fig. 8. It is known that for thermally stable crystal, C_P > C_V. As noted in the discussion for , for closed packed structure, C_V remains below the classical DP law. These two points thus suggest that the discrepancy seen in C_P at higher temperature can be attributed to γ_th and to some extent electronic part of C_V, as thermal expansion is already in confirmation with experimental results, cf. Eq. (11). Improving will also improve C_P. However, we believe that the major source of discrepancy in C_P is γ_th. Furthermore, extensive studies based on phonon frequency shifts for thermal anharmonicity in elemental metals^[42,50–52] all indicate that the largest anharmonicity in the crystal is related to thermal expansion. However, the present study contradicts this assertion near the melting temperature, and proposes non-negligible contribution from other than volume dilation mechanism.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. (color online) Variation of specific heat at constant pressure. Present result: lines (blue: An, red: QHA). Experimental results (symbol) are from Ref. [48].

Shown in Figs. 9 and 10 are the results for relative zero-pressure enthalpy or internal energy (H) and entropy difference ΔS = S(T)–S(0). Good agreement is observed for both these caloric properties with experimental data when anharmonicity is included. Thus, this observation confirms the importance of anharmonicity at higher temperature and low pressure. The discrepancy in enthalpy at high-T represents the sensitivity of this caloric property to the fitting parameters m and . We believe that it may be possible to determine m and by tuning both thermal expansion coefficient and enthalpy. Nonetheless, entropy seems to be less sensitive to specific choice of and m.

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. (color online) Relative enthalpy as a function of temperature. Present result: lines (blue: An, red: QHA). Experimental results (symbol) are from Ref. [48].

To address the discrepancies in enthalpy and C_P, we argue as follows. We propose that the thermal expansion is only insufficient to account for deducing correct caloric properties, e.g. C_P and enthalpy. This inconsistency requires systematic investigation of thermal response of crystal, which can be examined, for instance, in terms of Grüneisen parameter and thermal pressure. Figure 11 shows results for thermodynamic Grüneisen parameter. These results increase with temperature as melting is approached. Since it is difficult to hold the sample at constant volume experimentally at high-T, it puts restriction on the determination of explicit temperature variation of γ_th ≡ γ_th(T). Though all theoretical studies^[53] that compute γ_th(T) at constant volume suggest very weak temperature dependence, in practice γ_th(V₀(T)) evaluated at zero constant pressure condition increases with T due to combined effect of thermal dilation and its explicit temperature dependence.

Fig. 10.

	Figure Option ViewDownloadNew Window
	Fig. 10. (color online) Relative Entropy versus temperature. Present result: lines (blue: An, red: QHA). Experimental results (symbol) are from Ref. [48].

Fig. 11.

	Figure Option ViewDownloadNew Window
	Fig. 11. (color online) Temperature variation of thermodynamic Grüneisen parameter. Present results: lines (blue: An, red: QHA). Experimental results (star) and theoretical results (solid square) are quoted from Ref. [53].

To focus further on this anharmonic nature of the system, we show the product α_VB_T versus T using their equilibrium values in Fig. 12. By definition, the term α_VB_T denotes the thermal pressure P_th. Clearly, the graph in Fig. 12 can be divided into two temperature regions: the low temperature region in which P_th exhibits steep linear increase and the high temperature region. The slope in Fig. 12 alters its sign to (moderate) negative value, showing that this physical quantity decreases with temperature. We identify this inflection point also as the characteristic point (T_C) which again separates low-T quantum regime to high-T classical domain. It can be inferred from the graph of C_V versus T that at T_C, specific heat attains its classical value. Observation of Fig. 12 clearly suggests that up to T = T_C ≈ 200 K, thermal pressure P_th is positive, indicating that the absorbed thermal energy is responsible for rapid thermal expansion. It is important to note that T_C so determined is close to conventional Debye temperature and also to a temperature where ionic part of free energy is zero, shown in Fig. 3. In the language of thermodynamics, thermal energy does work on a system and results in an expansion below T_C. Further, since T_C so determined is usually low (about 1/3T_M or even lower) for majority class of materials, thermo-elasticity (e.g. Poisson ratio) is not severely affected by temperature below and slightly above T_C. In this situation, it is reasonable to assume that in the first approximation, the Poisson ratio is independent of temperature. It is natural to think that lattice vibrations are harmonic or near-so, such that thermal expansion can be well treated within the quasiharmonic phonon dynamics. However, at elevated temperature one expects departure from QH behaviour. For instance, at higher temperatures, increase in α_V is rather slow (Fig. 4), and thermal energy is absorbed to perform uncorrelated, faster, and large-amplitude oscillations, eventually resulting into varieties of effects like creation of vacancies, plasticity, dislocation, and wave-vector dependent softening of phonon-frequencies due to combined effect of volume expansion and temperature. Experimentally, these can be seen as line broadening in the observed spectra, more-so when sample is close to the melting temperature. As a plausible explanation, these intrinsic phenomena can be related to higher order potential derivatives^[49] and phonon scattering kinematics.^[54] Identifying the QHA as the first-order anharmonicity, the present anharmonic contribution F_An(V,T) is thought of as the second-order term. This is found to be insufficient to determine accurate caloric properties like enthalpy and C_P close to the melting temperature even when thermal expansion is accurately determined.

Fig. 12.

	Figure Option ViewDownloadNew Window
	Fig. 12. (color online) Temperature dependence of the product of volume thermal expansion coefficient and isothermal bulk modulus.

We report thermodynamic analysis for fcc silver up to melting temperature. The systematic evaluation of several thermal properties with and without including anharmonic contribution to ab initio phonon dynamics reveals the following important points. The electronic free energy is important only at high enough temperature, and vacancy contribution is ignored. Intrinsic phonon anharmonicity starts dominating above characteristic temperature T_C. The present study also reveals that F_An(V,T) ∼ T² correction is sufficient for correct evaluation of , P_th, and γ_th. Mutual comparison of for copper and silver suggests larger anharmonicity in silver. It is known that the decreases in the product α_VB_T and C_V (T) at high temperature suggest proper account of anharmonicity,^[8] and careful analysis of C_P and γ_th is carried out to elucidate the underlying mechanism. The discrepancy in caloric properties is attributed to higher order anharmonicity, which is important close to melting temperature. It can be related to higher order potential derivatives,^[49] an evaluation of which is intrinsically complicated from the first principles. This may include the contribution from vacancy formation. One possible way to tackle this pre-melting scenario is by keeping higher order perturbation terms for anharmonic free energy and treating their coefficients as fitting parameters. We conclude that the uncorrelated and large-amplitude lattice vibrations at high temperature raise dominating intrinsic thermal stress mechanism that surpasses the phonon-anharmonism and maybe vacancy contribution, which requires further attention.